#include <opencv2/opencv.hpp>
using namespace cv;

//主函数
int main()
{
    //①以灰度图像读取
    Mat srcimage = imread("/Users/lukeskywalker/Desktop/WechatIMG1.jpeg", 0);
    if (!srcimage.data)
    {
        printf("falsch!\n");
    }
    imshow("原始图像", srcimage);

    //2,将输入尺寸扩大到最佳尺寸，边界用0填充
    int m = getOptimalDFTSize(srcimage.rows);
    int n = getOptimalDFTSize(srcimage.cols);
    Mat padded;
    copyMakeBorder(srcimage, padded, 0, m - srcimage.rows, 0, n - srcimage.cols, BORDER_CONSTANT, Scalar::all(0)); //初始化元素为0，扩充

    //3,傅里叶变化结果（实部虚部）分配储存空间
    //将plannes数组组合合并成一个多通道的数组complexI
    Mat plannes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)}; //新建一个两页的array，其中第一页用扩展后的图像初始化，第二页初始化为0
    Mat complexI;
    merge(plannes, 2, complexI);

    //4，进行就地傅里叶变换
    dft(complexI, complexI);

    //5，将复数转化为幅值，即log（1+sqrt（Re(DFT(I)^2+Im(DFT(I))）^2）
    split(complexI, plannes); //多通道分离为几个单通道数组,plannes[0]=Re(DFT(I),plannes[1]=Im(DFT(I)
    magnitude(plannes[0], plannes[1], plannes[0]);
    Mat magnitudeImage = plannes[0];

    //6，进行对数尺度（logarithmic scale）缩放
    magnitudeImage += Scalar::all(1);    //等效于  magnitudeImage = magnitudeImage + Scalar::all(1);
    log(magnitudeImage, magnitudeImage); //求自然对数

    //7，剪切和重新分布幅度图象限
    //若有奇数行/列，进行频谱裁剪
    magnitudeImage = magnitudeImage(Rect(0, 0, magnitudeImage.cols & -2, magnitudeImage.rows & -2));
    //重新排列傅里叶图像的象限，是的原点位于图像中心
    int cx = magnitudeImage.cols / 2;
    int cy = magnitudeImage.rows / 2;
    Mat q0(magnitudeImage, Rect(0, 0, cx, cy)); //Roi区域的左上，右上，左下，右下
    Mat q2(magnitudeImage, Rect(0, cy, cx, cy));
    Mat q3(magnitudeImage, Rect(cx, cy, cx, cy));
    Mat q1(magnitudeImage, Rect(cx, 0, cx, cy));

    //交换象限（左上和右下）
    Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    //交换象限（右上与左下）
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);

    //8，归一化，用0到1间的浮点数将矩阵变化为可视的图像格式
    normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);

    //9,显示效果原图
    imshow("频谱幅值", magnitudeImage);
    waitKey();

    return 0;
}
